require(Rcpp)

today <- '20240725'
setwd(wkdir)
output_version <- Sys.Date()

##
## NTL AS SIZE VAR
##
ntl_list_fn <- list.files(ntldir, pattern='.tif',  full.names = FALSE)

#Tracks, improved, paved, highways and speed (r) but do Euclidean distance only (d) => two measures will be different if two measures are not strongly correlated. 
#Highway = 80. paved = 60. improved = 40. dirt = 12. no road = 6
#Sigma = 1, 2, 3, 4 (1 2 3 4)
#For example, baseline = lcmacro4 for each cell-year
#There are four surfaces:
#8 = highway (highways are paved roads but with at least 2 x 3 lanes) 
#1 = paved road
#3 = improved road
#5/0 = earthen road


## speeds in kph
roadcat2speed <- function(road_cat) {
  
  if (road_cat==8) {
    return(80)
  }
  if (road_cat==1) {
    return(60)
  }
  if (road_cat==3) {
    return(40)
  }
  if (road_cat==5) {
    return(12)
  } 
  if (road_cat==0) {
    return(6)
  } else {
    return(6)
  }
}

adm0_lake_chad_shp = paste0(wkdir,'/local/gis_data/003_boundaries/lake_chad/lake_chad_adm0')
adm0_reg_shp <- rgdal::readOGR(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
adm0_lake_chad_sf <- read_sf(dsn=dirname(adm0_lake_chad_shp), layer=basename(adm0_lake_chad_shp))
fishnet_ext <- extent(adm0_reg_shp)

##
## FISHNET
##
fishnet_lake_chad_ntl_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_x_d01')
fishnet_only <- read_sf(dirname(fishnet_lake_chad_ntl_x_d01_shp),basename(fishnet_lake_chad_ntl_x_d01_shp))
fishnet_only$CMR[is.na(fishnet_only$CMR)]<-0
fishnet_only$NER[is.na(fishnet_only$NER)]<-0
fishnet_only$NGA[is.na(fishnet_only$NGA)]<-0
fishnet_only$TCD[is.na(fishnet_only$TCD)]<-0

fishnet_tif <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_x_d01.tif')
fishnet_grid <- raster(fishnet_tif)

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
countries_GID_0 <- unique(ntl.zones$GID_0)
# add area km2 #
ntl.zones$sf_area_km2 <- st_area(ntl.zones) / 1000000

setwd(wkdir)

##
## roads
##
countries <-  c( 'Niger', 'Nigeria', 'Chad', 'Cameroon' )


# slope
Slope_globe <- raster(file.path(slope_dir,"meanslope.tif"))
Slope_crop <- crop(Slope_globe, alignExtent(fishnet_grid,Slope_globe,snap='out'))
extent(Slope_crop)
extent(fishnet_grid)
Slope_radians = Slope_crop / 100
Slope <- Slope_radians*pi/180
rm(Slope_globe)

##
## return ((6 * exp(-3.5 * abs(x + 0.05))) * 0.6)
## waldo in km/h
##
offpath <- function(x) {
  return ((6 * exp(-3.5 * abs(x*pi/180 + 0.05))) * 0.6)
}

## radians to degrees
onpath <- function(x) {
  return (6 * exp(-3.5 * abs(tan(x) + 0.05) ) )
}

grid_res_km = 11.1

# input: slope in degrees 
# output: km / h 
slope_kph <- raster::calc(Slope, onpath)

# cost in hours per cell
# time = distance / speed
# raster::area no weights
slope_hr.prj <- projectRaster(slope_kph,fishnet_grid, res = 0.1, method='bilinear')
slope_hr.prj <- mask(slope_hr.prj,fishnet_grid)
#slope_hr.cell = sqrt(raster::area(slope_hr.prj)) / slope_hr.prj

require(lwgeom)
road_lk_chad_shp = paste0(wkdir,'/local/gis_data/018_transportation/roads_panel')

if (!file.exists(paste0(road_lk_chad_shp,'.shp'))){
  roadsafr_sf <- st_read(dsn=roadsafr_dir,layer='rds_remi_afr') %>% 
    sf::st_make_valid()
  
  reg_sf <- st_read(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
  GID_0 <- reg_sf$GID_0
  
  roadsafr_prj_sf <- st_transform(roadsafr_sf,crs(reg_sf))
  roads <- subset(roadsafr_prj_sf, COUNTRY %in% countries)
  
  write_sf(roads,dirname(road_lk_chad_shp),basename(road_lk_chad_shp),driver='ESRI Shapefile')
} else {
  roads <- read_sf(dirname(road_lk_chad_shp),layer=basename(road_lk_chad_shp))
  roads = roads %>% dplyr::mutate(dTCD = ifelse(COUNTRY == "Chad", 1,0))
  roads = roads %>% dplyr::mutate(dCMR = ifelse(COUNTRY == "Cameroon", 1,0)) 
  roads = roads %>% dplyr::mutate(dNER = ifelse(COUNTRY == "Niger", 1,0)) 
  roads = roads %>% dplyr::mutate(dNGA = ifelse(COUNTRY == "Nigeria", 1,0))
  roads = roads %>% dplyr::mutate(ISO3 = ifelse(COUNTRY == "Nigeria", "NGA",
                                                ifelse(COUNTRY == "Niger", "NER",
                                                              ifelse(COUNTRY == "Cameroon", "CMR",
                                                                     ifelse(COUNTRY == "Chad", "TCD",NA)))))
}   

##
## timecost ALL PIXELS
##
timecost_t <- function(ryear_t,iso3roadlist) {
  # include only countries of interest
  # roads
  roads.t = subset(roads,ISO3 %in% iso3roadlist)
  iso3roads <- str_c(sort(unique(roads.t$ISO3)), collapse = "")
  
  # by period
  roads.t = subset(roads.t,select=c(COUNTRY,eval(as.name(paste0("R",ryear_t)))))
  
  fishnet_grid_c_fn = paste0(wkdir,'/local/gis_data/003_boundaries/fishnet_',iso3roads,'.tif')
  if(!file.exists(fishnet_grid_c_fn)){
    fishnet_c <- subset(fishnet_only,select=c("ID",iso3roadlist))
    cond = str_c(paste0(iso3roadlist," == 1"),collapse=" | ")
    fishnet_c <- fishnet_c %>% filter(eval(parse(text = cond)))
    fishnet_c_sp <-as(st_as_sf(fishnet_c),Class="Spatial")
    fishnet_grid_c <- mask(fishnet_grid,fishnet_c_sp)
    writeRaster(fishnet_grid_c, fishnet_grid_c_fn, driver= "GTiff",overwrite=TRUE)
  } else {
    fishnet_grid_c <- raster(fishnet_grid_c_fn)
  }
  # max speed for each cell 
  roads.t$speed = unlist(lapply(unlist(st_drop_geometry(roads.t[paste0("R",ryear_t)])),FUN=roadcat2speed))
  require(stars)
  roads.grd.fn = paste0(wkdir,'/local/gis_data/018_transportation/road_max',iso3roads,ryear_t,'.tif')
  if(!file.exists(roads.grd.fn)){
    roads.grd.t = raster::rasterize(as(roads.t,Class="Spatial"),fishnet_grid,field = 'speed', fun='max', background=NA)
    writeRaster(roads.grd.t, roads.grd.fn, driver= "GTiff",overwrite=TRUE)
  } else {
    roads.grd.t <- raster(roads.grd.fn)
  }
  
  roads.hr =  roads.grd.t  
  # slope_hr.prj kph in cell
  slope_hr.mask <- mask(slope_hr.prj,fishnet_grid_c)
  costgrid.stack <- stack(slope_hr.mask,roads.hr)
  # maximum speed of road path and natural path in a cell
  cost.grd.t <- stackApply(costgrid.stack, c(1,1), fun='max',na.rm=TRUE) 
  
  # mask country only including full border cell #
  cost.grd.t <- mask(cost.grd.t,fishnet_grid_c)
  
  timecost <- paste0(wkdir,'/local/gis_data/018_transportation/timecostgrid_min',iso3roads,ryear_t,'.tif')
  if (!file.exists(timecost)) {
    writeRaster(cost.grd.t,timecost,driver='GTiff',overwrite=TRUE)
  }
  
  ## update cost-distance grid with background levels from waldo tobler function
  ## the values in the input raster are costs. 
  ## take the arithmetic mean of the cost first and then take the reciprocal of that to get the conductance. 
  ## ie. The traveller experiences half of the cost of the origin cell and half of the cost of the destination cell, (cost1 + cost2)/2.
  return(cost.grd.t)
}


##
## timecost NO BH PIXELS
## 
timecost_BH_t <- function(ryear_t,iso3roadlist,bhbr_kph) {
  roads.t = subset(roads,ISO3 %in% iso3roadlist)
  iso3roads <- str_c(sort(unique(roads.t$ISO3)), collapse = "")
  
  # by period
  roads.t = subset(roads.t,select=c(COUNTRY,eval(as.name(paste0("R",ryear_t)))))
  
  fishnet_adm0_only$d_road <- lengths(st_intersects(fishnet_adm0_only,roads.t))
  fishnet_roads_out <- fishnet_adm0_only %>%
    st_drop_geometry() %>%
    dplyr::select(OBJECTID,d_road) %>%
    dplyr::mutate(d_road = ifelse(d_road>0,1,0),
                  year = ryear_t) %>%
    rename_with(tolower, .cols = everything()) %>%
    as.data.frame()
  fishnet_roads_out_fn <- file.path(wkdir,'proc_data',sprintf('ROADS%syn.dta',ryear_t))
  haven::write_dta(fishnet_roads_out,fishnet_roads_out_fn)
  
  # BH area 5kph only
  fishnet_grid_c_fn = paste0(wkdir,'/local/gis_data/003_boundaries/fishnet_',iso3roads,'.tif')
  if(!file.exists(fishnet_grid_c_fn)){
    fishnet_c <- subset(fishnet_only,select=c("ID",iso3roadlist))
    cond = str_c(paste0(iso3roadlist," == 1"),collapse=" | ")
    fishnet_c <- fishnet_c %>% filter(eval(parse(text = cond)))
    fishnet_c_sp <-as(st_as_sf(fishnet_c),Class="Spatial")
    fishnet_grid_c <- mask(fishnet_grid,fishnet_c_sp)
    writeRaster(fishnet_grid_c, fishnet_grid_c_fn, driver= "GTiff",overwrite=TRUE)
  } else {
    fishnet_grid_c <- raster(fishnet_grid_c_fn)
  }
  
  # transform road classes into travel cost to cross a pixel
  # grid resolution in km
  
  # max speed for each cell 
  roads.t$speed = unlist(lapply(unlist(st_drop_geometry(roads.t[paste0("R",ryear_t)])),FUN=roadcat2speed))
  summary(roads.t$speed)
  require(stars)
  roads.grd.fn = paste0(wkdir,sprintf('/local/gis_data/018_transportation/road_max_bh%skph',bhbr_kph),iso3roads,ryear_t,'.tif')
  if(!file.exists(roads.grd.fn)){
    # intersection of roads in BH area
    bhbr_lines <- bhbr %>%
      dplyr::select(geometry) %>%
      mutate(BHBR = 1) %>%
      st_cast(.,"LINESTRING")
    intersections <- st_intersection(roads.t, bhbr_lines) %>% st_cast(.,"POINT")
    roads.t_split <- st_split(roads.t, intersections) %>%
      st_collection_extract(.,"LINESTRING")
    
    # assign lines inside BHBR
    roads.t_split$within_bhbr <- as.numeric(st_within(roads.t_split, bhbr))
    roads.t_split$within_bhbr <- as.numeric(ifelse(is.na(roads.t_split$within_bhbr), 0, 1))
    table(roads.t_split$within_bhbr)
    # roads in BHBR are set to lowest speed at 6kph
    roads.t_split$speed[roads.t_split$within_bhbr==1] <- bhbr_kph
    roads.grd.t = raster::rasterize(as(roads.t_split,Class="Spatial"),fishnet_grid,field = 'speed', fun='max', background=NA)
    ## modify values in BH area to background speed 6kph
    writeRaster(roads.grd.t, roads.grd.fn, driver= "GTiff",overwrite=TRUE)
  } else {
    roads.grd.t <- raster(roads.grd.fn)
  }
  
  roads.hr =  roads.grd.t
  slope_hr.mask <- mask(slope_hr.prj,fishnet_grid_c)
  costgrid.stack <- stack(slope_hr.mask,roads.hr)
  
  # least cost of road path and natural path
  cost.grd.t <- stackApply(costgrid.stack, c(1,1), fun='max',na.rm=TRUE) 
  
  # mask country only including full border cell #
  cost.grd.t <- mask(cost.grd.t,fishnet_grid_c)
  
  timecost <- paste0(wkdir,sprintf('/local/gis_data/018_transportation/timecostgrid_min_bhbr%skph',bhbr_kph),iso3roads,ryear_t,'.tif')
  if (!file.exists(timecost)) {
  writeRaster(cost.grd.t,timecost,driver='GTiff',overwrite=TRUE)
  }
  
  ## update cost-distance grid with background levels from waldo tobler function
  ## the values in the input raster are costs. 
  ## take the arithmetic mean of the cost first and then take the reciprocal of that to get the conductance. 
  return(cost.grd.t)
}


##
## END COST
##
s <- stack(fishnet_grid)
iso3list <- c("CMR","NGA","NER","TCD")

#Tracks, improved, paved, highways and speed (r) but do Euclidean distance only (d) => two measures will be different if two measures are not strongly correlated. 
#Highway = 80. paved = 60. improved = 40. dirt = 12. no road = 6
#Sigma = 1, 2, 3, 4 (1 2 3 4)

calcMarketAccessParallel <<- function(cost.c.t,fromCoordsWeight,sizentlgrid,trade_param,excl_km,.packages=c("gdistance","raster","igraph", "Matrix")){
  
  sizentlgrid = extend(sizentlgrid,cost.c.t,value=NA)
  
  ncoords<-dim(fromCoordsWeight)[1]
  
  cluster <- makeCluster(18, type = "SOCK")
  registerDoSNOW(cluster)
  
  ## parallel access of conflict calculation with trade parameter alpha ##
  ## for each conflict event location in grid cell by time-step ##
  m<-foreach(n = 1:ncoords, 
             .packages=c("gdistance","raster","igraph", "Matrix") ) %dopar% {
               require(gdistance)
               fromCoordsMat <- c(fromCoordsWeight[n,1],fromCoordsWeight[n,2])
               acc_cost.hour <- gdistance::accCost(cost.c.t, fromCoordsMat)
               fun1 <- function(x){ifelse((x < 1), 1, x)}
               acc_cost.hour[acc_cost.hour!=0] <- fun1(acc_cost.hour[acc_cost.hour!=0])
               excldistkm = raster::distanceFromPoints(acc_cost.hour, fromCoordsMat)
               fun <- function(x) {ifelse((x <= (excl_km * 1000)), 0, 1)}
               excl20km = raster::calc(excldistkm,fun)
               size_cell_i = as.double(fromCoordsWeight[n,3])
               # exclude pixel center within 20km of center point #
               acc_cost.hour[acc_cost.hour!=0] <- (size_cell_i)  * ((excl20km[acc_cost.hour!=0]) * (acc_cost.hour[acc_cost.hour!=0]^-trade_param))
               return(acc_cost.hour)
             }
  stopCluster(cluster)
  
  # sum over all points
  cum_stack <- stack(m)
  # add 1 so it's possible to take ln #
  cum_sum <- sum(cum_stack)
  cum_sum_mask <- mask(cum_sum,sizentlgrid)
  return(cum_sum_mask)
}

ntl_year=c(2008)
rd_years = 2008


##
## Travel time by region
## Regional level open borders
##
## Market access measures to other NTL (different parameters): 
## (i) in same country only, 
## (ii) in all countries not driving through BH area
##

##
## tt with no bhbr area
##
##
## BH between rivers
##
bhbr_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/bhbr_dst')
bhbr <- read_sf(dsn=dirname(bhbr_shp),layer=basename(bhbr_shp))

## rasterize uses center cell
bhbr_grid <- rasterize(bhbr,fishnet_grid,field="OBJECTID",fun=mean)
bhbr_grid[!is.na(bhbr_grid)] <- -1
bhbr_grid[is.na(bhbr_grid)] <- 1


##
## o=open border
##
border_sf = st_read(dsn=paste0(wkdir,'/local/gis_data/003_boundaries'),layer='gadm36_lake_chad_ctr_0_ln')
border_sf["borderid"]=1
border_fn=paste0(wkdir,'/local/gis_data/018_transportation/road_borders_fishnet.tif')
if(!file.exists(border_fn)){
  border <- rasterize(border_sf,fishnet_grid,field="borderid",fun=mean,background=0)
  writeRaster(border,border_fn)
} else{
  border <- raster(border_fn)
}

## need border points snapped to roads - manual
bpost_sf <- st_read(paste0(wkdir,'/local/gis_data/003_boundaries/border_post.shp'))
bpost_sf["bid"]=1
bpost <- rasterize(bpost_sf,fishnet_grid,field="bid",fun=mean,background=0)

# add border and border post dummy 
bopen = border + bpost
bopen[bopen==1]=NA

# copy bopen for 4 hour wait
bopen4 = bopen
bopen[bopen==2]=1
bopen[bopen==0]=1

# add border cost time in hours
bopen4[bopen4==2] = 4

ntl_years <- c(2008:2008)

outdir <-'D:/temp'

ntl_years <- c(2008:2008)

## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

for (ntlyear_t in ntl_years){
  out_fn <- paste0(wkdir,'/proc_data/ma_ntl_excl20km_bhbr_shocks',ntlyear_t,today,'.dta')
  #if(!file.exists(out_fn)) {
      s <- stack(fishnet_grid)
      
      ryear_t = 2008
      print(paste0("calc regional MA with NTL: ",ntlyear_t," roads ",ryear_t))  
      
      ##
      ## NEW
      ##
      
      print("porous border loses significane")
      cost.grd.4B_P_bhbr6kph = timecost_BH_t(ryear_t,c("NER","TCD","CMR","NGA"),6)
      cost.grd.4B_P_bhbr0kph = timecost_BH_t(ryear_t,c("NER","TCD","CMR","NGA"),0)
      cost.grd.4B_P_bhbr_all = timecost_t(ryear_t,c("NER","TCD","CMR","NGA"))
      
      print("added open border + 4 hours")
      cost.grd.4B_O4_bhbr6kph <- cost.grd.4B_P_bhbr6kph + bopen4
      cost.grd.4B_O4_bhbr0kph <- cost.grd.4B_P_bhbr0kph + bopen4
      cost.grd.4B_O4_bhbr_all <- cost.grd.4B_P_bhbr_all + bopen4
           
      bhlist = c("ma_4S6_P","ma_4S0_P","ma_4S9_P",
                 "ma_4S6_O4","ma_4S0_O4","ma_4S9_O4")
      
      for (bh_i in bhlist){
          print(sprintf("Roads scenario %s", bh_i))
          
           bordervar = strsplit(bh_i, "_")[[1]][3]
           cval = strsplit(bh_i, "_")[[1]][2]
          if (bh_i=="ma_4S6_P"){ cost.grd.reg.t = cost.grd.4B_P_bhbr6kph
              cvar="CMRNERNGATCD"}
          if (bh_i=="ma_4S0_P"){ cost.grd.reg.t = cost.grd.4B_P_bhbr0kph
              cvar="CMRNERNGATCD"}
          if (bh_i=="ma_4S9_P"){ cost.grd.reg.t = cost.grd.4B_P_bhbr_all
            cvar="CMRNERNGATCD"}
          
           # open border + 4 hours wait
          if (bh_i=="ma_4S6_O4"){ cost.grd.reg.t = cost.grd.4B_O4_bhbr6kph
           cvar="CMRNERNGATCD"}
          if (bh_i=="ma_4S0_O4"){ cost.grd.reg.t = cost.grd.4B_O4_bhbr0kph
           cvar="CMRNERNGATCD"}
          if (bh_i=="ma_4S9_O4"){ cost.grd.reg.t = cost.grd.4B_O4_bhbr_all
           cvar="CMRNERNGATCD"}
           
           fishnet_grid_c_fn = paste0(wkdir,'/local/gis_data/003_boundaries/fishnet_',cvar,'.tif')
           fishnet_grid_c <- raster(fishnet_grid_c_fn)
           
           ## MA in country
           print('NTL size var')
           # return 1st occurrence of NTL in year #
           ntl.year.stack.index <- min(which(as.integer(substr(ntl_list_fn,4,7)) %in% ntlyear_t))
           ntl.year <- raster(paste(ntldir, ntl_list_fn[ntl.year.stack.index],sep='/'))
           
           ## crop to country
           ntl_c_t <- crop(ntl.year,fishnet_grid_c,snap='near')
           
           ## aggregate points to master fishnet
           ntlfactor = 12 
           
           ## SUM OF NTL
           ntlval_c_t_agg <- aggregate(ntl_c_t,fact=ntlfactor,fun=sum,background=0)
           
           ## overlay of nighttime lights directly
           fishnet_out_fn <- file.path(wkdir,'proc_data',sprintf('NTLSUM%syn.dta',ntlyear_t))
           if(!file.exists(fishnet_out_fn)){
             fishnet_adm0_only$ntlsumagg <- exact_extract(ntlval_c_t_agg, fishnet_adm0_only, 'sum')
             
             #convert to dummy of NTL in cell or NOT
             fishnet_out <- fishnet_adm0_only %>%
               st_drop_geometry() %>%
               dplyr::select(OBJECTID,ntlsumagg) %>%
               dplyr::mutate(d_ntlsumagg = ifelse(ntlsumagg>0,1,0),
                             year = ntlyear_t) %>%
               rename_with(tolower, .cols = everything())
             
             haven::write_dta(fishnet_out,fishnet_out_fn)
           }
           
          # convert from km/hr cost to m/hr cost (implicit transition correction units)
          cost.grd.reg.t = mask(cost.grd.reg.t,fishnet_grid_c) * 1000
          print("cost.grd.reg.t = mask(cost.grd.reg.t,fishnet_grid_c) * 1000 / 60")
          
          # not permeability
          cost.c.t = gdistance::transition(cost.grd.reg.t, transitionFunction=function(x) mean(x), directions=8)
          cost.c.t <- geoCorrection(cost.c.t, type="c")

          ##
          ## add variation of NTL = 0 in BHBR or not
          ##
          ntl_bhbr <- c('bh0','bh1')
          
          for(ntl_bhbr_i in ntl_bhbr){
            
            # start fresh with NTL value
            # mask to region (1C,3C,4C,5C)
            ntlval_c_t <- mask(ntlval_c_t_agg,fishnet_grid_c)  
            
            print(ntl_bhbr_i)
            
            ## only certain scenarios ##
            if((bh_i=="ma_4S9_P" & ntl_bhbr_i=="bh1") | (bh_i=="ma_4S9_O4" & ntl_bhbr_i=="bh1")){
              print(sprintf("DO NOT skip iteration of %s for %s",ntl_bhbr_i,bh_i))
              print("BASELINE SCENARIO NO SHOCKS FOR COMPARISON with run in 2021")
              #next()
            }
            if((bh_i=="ma_4S0_P" & ntl_bhbr_i=="bh1") | (bh_i=="ma_4S0_O4" & ntl_bhbr_i=="bh1")) {
              print(sprintf("skip iteration of %s for %s",ntl_bhbr_i,bh_i))
              next()
            } 
            if((bh_i=="ma_4S6_P" & ntl_bhbr_i=="bh0") | (bh_i=="ma_4S6_O4" & ntl_bhbr_i=="bh0")){
              print(sprintf("skip iteration of %s for %s",ntl_bhbr_i,bh_i))
              next()
            } 
            
            
            if(ntl_bhbr_i=="bh0"){
              # bhbr region = 0 ntl
              print("setting NTL values to 0 in BHBR area")
              bhbr_grid0 <- bhbr_grid
              # grid is 0 in BHBR area and 1 elsewhere
              bhbr_grid0[bhbr_grid0==-1]=0
              
              # edit ntl values to 0
              ntlval_c_t <- ntlval_c_t * bhbr_grid0 
              
            }
            
            
            if (cellStats(ntlval_c_t,stat = 'max')!=0){
              
              ## sizevar is NTL
              fromcoordsNTLWeight = rasterToPoints(ntlval_c_t,values=TRUE)
              fromcoordsNTLWeight = fromcoordsNTLWeight[fromcoordsNTLWeight[,3]>0,]
              
              ##
              ## calculate excluding 0
              ##
              excl_list <- c(0)
              for (exclvar in excl_list) {
                alist=c(1,2,3,4) 
                for (alpha in alist){
                  print(paste0("BHBR NTL vals: ",ntl_bhbr_i," alpha : ",alpha, " and excluding buffer (km) :", exclvar))
                  acc_ma_t <- calcMarketAccessParallel(cost.c.t,fromcoordsNTLWeight,ntlval_c_t,alpha,exclvar)
                  acc_ma_t <- extend(acc_ma_t,s[[1]])
                  ## [iso3] ma=market access; r=roads,c=closed border,alpha,ntl year ##
                  ntlyear_2d_t <- sprintf('%02d', ntlyear_t %% 100)
                  decayvar = "tt"
                  names(acc_ma_t) <- paste("unlogntlma",ntlyear_2d_t,cval,decayvar,bordervar,sprintf('%02d', exclvar),alpha,ntl_bhbr_i,sep="_")
                  s<-stack(s,acc_ma_t)
                }
              }
            }
            
          }
          
          
      }
  
  # convert stack to DF
  out_df <- as.data.frame(s)
  out_df <- out_df %>% dplyr::select(sort(names(.)))
  names(out_df)
  attr(out_df$fishnet_lake_chad_ntl_x_d01, "label") <- "gridid01 identifier"
  
  print("bh0 is ntl values in BHBR are set to 0")
  attr(out_df$unlogntlma_08_4S0_tt_P_00_1_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR nat. excl 0km BHBR NTL 0 a=1"
  attr(out_df$unlogntlma_08_4S0_tt_P_00_2_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR nat. excl 0km BHBR NTL 0 a=2"
  attr(out_df$unlogntlma_08_4S0_tt_P_00_3_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR nat. excl 0km BHBR NTL 0 a=3"
  attr(out_df$unlogntlma_08_4S0_tt_P_00_4_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR nat. excl 0km BHBR NTL 0 a=4"
  attr(out_df$unlogntlma_08_4S0_tt_O4_00_1_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR nat. excl 0km BHBR NTL 0 a=1"
  attr(out_df$unlogntlma_08_4S0_tt_O4_00_2_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR nat. excl 0km BHBR NTL 0 a=2"
  attr(out_df$unlogntlma_08_4S0_tt_O4_00_3_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR nat. excl 0km BHBR NTL 0 a=3"
  attr(out_df$unlogntlma_08_4S0_tt_O4_00_4_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR nat. excl 0km BHBR NTL 0 a=4"
  
    
  ##
  ##
  print("bh0 is ntl values in BHBR are set to 0")
  attr(out_df$unlogntlma_08_4S9_tt_P_00_1_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR all excl 0km BHBR NTL 0 a=1"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_2_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR all excl 0km BHBR NTL 0 a=2"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_3_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR all excl 0km BHBR NTL 0 a=3"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_4_bh0, "label") <- "unlog ntl MA 2008 Porous border BHBR all excl 0km BHBR NTL 0 a=4" 
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_1_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR all excl 0km BHBR NTL 0 a=1"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_2_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR all excl 0km BHBR NTL 0 a=2"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_3_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR all excl 0km BHBR NTL 0 a=3"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_4_bh0, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR all excl 0km BHBR NTL 0 a=4"
  
  print("bh1 is ntl values in BHBR are as-is and bad roads")
  attr(out_df$unlogntlma_08_4S6_tt_P_00_1_bh1, "label") <- "unlog ntl MA 2008 BHBR 6kph excl  0km BHBR NTL ALL a=1"
  attr(out_df$unlogntlma_08_4S6_tt_P_00_2_bh1, "label") <- "unlog ntl MA 2008 BHBR 6kph excl  0km BHBR NTL ALL a=2"
  attr(out_df$unlogntlma_08_4S6_tt_P_00_3_bh1, "label") <- "unlog ntl MA 2008 BHBR 6kph excl  0km BHBR NTL ALL a=3"
  attr(out_df$unlogntlma_08_4S6_tt_P_00_4_bh1, "label") <- "unlog ntl MA 2008 BHBR 6kph excl  0km BHBR NTL ALL a=4"
  

  attr(out_df$unlogntlma_08_4S6_tt_O4_00_1_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR 6kph excl 0km BHBR NTL ALL a=1"
  attr(out_df$unlogntlma_08_4S6_tt_O4_00_2_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR 6kph excl 0km BHBR NTL ALL a=2"
  attr(out_df$unlogntlma_08_4S6_tt_O4_00_3_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR 6kph excl 0km BHBR NTL ALL a=3"
  attr(out_df$unlogntlma_08_4S6_tt_O4_00_4_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR 6kph excl 0km BHBR NTL ALL a=4"
  
  print("bh1 has ntl values in BHBR are as-is")
  attr(out_df$unlogntlma_08_4S9_tt_P_00_1_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=1"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_2_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=2"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_3_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=3"
  attr(out_df$unlogntlma_08_4S9_tt_P_00_4_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=4"

  attr(out_df$unlogntlma_08_4S9_tt_O4_00_1_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR NO SHOCKS a=1"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_2_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR NO SHOCKS a=2"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_3_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR NO SHOCKS a=3"
  attr(out_df$unlogntlma_08_4S9_tt_O4_00_4_bh1, "label") <- "unlog ntl MA 2008 Open border +4hr BHBR NO SHOCKS a=4"
  
    
  # output market shocks file
  attr(out_df, "label") <- sprintf("prep_ma_ntl_tt_excl20km_server_ts_2024-07-23_wp.R last run on %s",Sys.Date())
  haven::write_dta(out_df,out_fn)  

  ##
  ##
  crosswalk_fn <- file.path(wkdir,'proc_data','CROSSWALKsorted.dta')
  crosswalk <- haven::read_dta(crosswalk_fn)
  
  out_obj <- merge(out_df,crosswalk,
                   by.x='fishnet_lake_chad_ntl_x_d01',
                   by.y='gridid01',
                   all.x=F)
  dim(out_obj)
  names(out_obj)
  out_obj <- out_obj %>%
    dplyr::select(!fishnet_lake_chad_ntl_x_d01) %>%
    relocate(objectid) %>%
    arrange(objectid)
  
  
  ##
  ## NEW LABELS 
  ##
  out_obj_fn <- paste0(wkdir,'/proc_data/ma_ntl_excl20km_bhbr_shocks_newlabels',ntlyear_t,today,'.dta')
  
  attr(out_obj$unlogntlma_08_4S0_tt_P_00_1_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 & BH RDS 0 a=1"
  attr(out_obj$unlogntlma_08_4S0_tt_P_00_2_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 & BH RDS 0 a=2"
  attr(out_obj$unlogntlma_08_4S0_tt_P_00_3_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 & BH RDS 0 a=3"
  attr(out_obj$unlogntlma_08_4S0_tt_P_00_4_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 & BH RDS 0 a=4"

  attr(out_obj$unlogntlma_08_4S0_tt_O4_00_1_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 & BH RDS 0 a=1"
  attr(out_obj$unlogntlma_08_4S0_tt_O4_00_2_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 & BH RDS 0 a=2"
  attr(out_obj$unlogntlma_08_4S0_tt_O4_00_3_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 & BH RDS 0 a=3"
  attr(out_obj$unlogntlma_08_4S0_tt_O4_00_4_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 & BH RDS 0 a=4"
  
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_1_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 a=1"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_2_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 a=2"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_3_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 a=3"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_4_bh0, "label") <- "unlog ntl MA 08 fix BH NTL 0 a=4"

  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_1_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 a=1"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_2_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 a=2"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_3_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 a=3"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_4_bh0, "label") <- "unlog ntl MA 08 open border +4hrs fix BH NTL 0 a=4"
  
  
  attr(out_obj$unlogntlma_08_4S6_tt_P_00_1_bh1, "label") <- "unlog ntl MA 08 fix BH RDS 0 a=1"
  attr(out_obj$unlogntlma_08_4S6_tt_P_00_2_bh1, "label") <- "unlog ntl MA 08 fix BH RDS 0 a=2"
  attr(out_obj$unlogntlma_08_4S6_tt_P_00_3_bh1, "label") <- "unlog ntl MA 08 fix BH RDS 0 a=3"
  attr(out_obj$unlogntlma_08_4S6_tt_P_00_4_bh1, "label") <- "unlog ntl MA 08 fix BH RDS 0 a=4"
  
  attr(out_obj$unlogntlma_08_4S6_tt_O4_00_1_bh1, "label") <- "unlog ntl MA 08 open border +4hrs fix BH RDS 0 a=1"
  attr(out_obj$unlogntlma_08_4S6_tt_O4_00_2_bh1, "label") <- "unlog ntl MA 08 open border +4hrs fix BH RDS 0 a=2"
  attr(out_obj$unlogntlma_08_4S6_tt_O4_00_3_bh1, "label") <- "unlog ntl MA 08 open border +4hrs fix BH RDS 0 a=3"
  attr(out_obj$unlogntlma_08_4S6_tt_O4_00_4_bh1, "label") <- "unlog ntl MA 08 open border +4hrs fix BH RDS 0 a=4"
  
  print("SCENARIO 10 IS TO DOUBLE CHECK VALUES WITH MA_NTL_TT_HR08.dta MODIFIED 2022-03-21")
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_1_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=1"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_2_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=2"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_3_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=3"
  attr(out_obj$unlogntlma_08_4S9_tt_P_00_4_bh1, "label") <- "unlog ntl MA 2008 BHBR NO SHOCKS a=4"

  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_1_bh1, "label") <- "unlog ntl MA 2008 Open border +4hrs BHBR NO SHOCKS a=1"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_2_bh1, "label") <- "unlog ntl MA 2008 Open border +4hrs BHBR NO SHOCKS a=2"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_3_bh1, "label") <- "unlog ntl MA 2008 Open border +4hrs BHBR NO SHOCKS a=3"
  attr(out_obj$unlogntlma_08_4S9_tt_O4_00_4_bh1, "label") <- "unlog ntl MA 2008 Open border +4hrs BHBR NO SHOCKS a=4"
  
  names(out_obj)
  
  haven::write_dta(out_obj,out_obj_fn)  
  
  ##
  ## produce dummy table with presence of NTL AND ROAD in year t
  ##
  fishnet_roads_out_fn <- file.path(wkdir,'proc_data',sprintf('ROADS%syn.dta',ryear_t))
  fishnet_d_roads <- haven::read_dta(fishnet_roads_out_fn)
  fishnet_out_fn <- file.path(wkdir,'proc_data',sprintf('NTLSUM%syn.dta',ntlyear_t))
  fishnet_d_ntl <- haven::read_dta(fishnet_out_fn)
  
  # merge
  fishnet_road_and_ntl <- merge(fishnet_d_roads,fishnet_d_ntl) %>%
    arrange(objectid,year)
  head(fishnet_road_and_ntl)
  fishnet_road_and_ntl <- fishnet_road_and_ntl %>%
    mutate(d_road_and_ntl = ifelse(d_road==1 & d_ntlsumagg==1,1,0)) %>%
    dplyr::select(objectid,d_road_and_ntl,year)
  
  attr(fishnet_road_and_ntl$d_road_and_ntl,"label") <- "==1 Roads and NTL present in cell"
  attr(fishnet_road_and_ntl,"label") <- sprintf("prep_ma_ntl_tt_excl20km_server_ts_2024-07-23_wp.R last run on %s",Sys.Date())
  fishnet_road_and_ntl_fn <- file.path(wkdir,'proc_data','NTLROADyn.dta')
  haven::write_dta(fishnet_road_and_ntl,fishnet_road_and_ntl_fn)
}

